**** Synthetic control placebo testing 

matrix deterrence = (0,0)

*** Load in data, loop over different controls
cd /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/botox/synth/100pct/
use growthgroupspanel.dta, clear 
drop if group == 999
levelsof group, local(placebos)

foreach p of local placebos{

	cd /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/botox/synth/100pct/placebos

	*** Grab full data and remove treated unit
	use /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/botox/synth/100pct/growthgroupspanel.dta, clear
	drop if group == 999 


	******************************************
	** Swap in placebos
	disp "Placebo group `p'"
	replace group = 999 if group == `p' 
	save growthgroupspanel_temp.dta, replace 

	*********************************************************************
	** Regress treated unit on lags to find best fit  

	levelsof group, local(levels)
	keep group dt pmt_amt 

	reshape wide pmt_amt, i(dt) j(group)
	tsset dt 


	matrix distanceLag = (0,0,0)
	matrix distanceLead = (0,0,0)
	* 3 columns: group, lag period, distance

	set matsize 11000


	*** Only want to keep botox pre-period
	local treatment_period = mofd(mdy(6,5,2007))
	replace pmt_amt999 = . if dt >= `treatment_period'

	** Find max lead so that there are 5 years of post-period data 
	summarize dt 
	local maxlead = r(max) - `treatment_period' + 1 - 60  
		** Max - treatment  + 1 = length of end window 
		** Take 60 off and that's the max you can shift before running out of data

	** Find max lag so that there are 3 years of pre-period data 
	summarize dt 
	local maxlag = `treatment_period' - r(min) - 36
		** Treatment period - min  = length of end window 


	gen diff = 0 
	*** Regress, storing leads and lags
	foreach g of local levels {
		disp `g'
		* Skip Treated group
		if `g' == 999{
				continue
			}
		foreach t of numlist 0/`maxlag'{
			*disp `t'
			qui replace diff = (pmt_amt999 - L`t'.pmt_amt`g')^2
			qui summarize diff 
			matrix distanceLag = (distanceLag \ `g', `t', r(sum)/r(N))
			}

		foreach t of numlist 0/`maxlead'{
			* disp `t'
			qui replace diff = (pmt_amt999 - F`t'.pmt_amt`g')^2
			qui summarize diff 
			matrix distanceLead = (distanceLead \ `g', `t', r(sum)/r(N))
			}
	}

	svmat distanceLag
	keep distanceLag*
	svmat distanceLead

	rename distanceLag1 laggroup
	rename distanceLag2 lagperiod
	rename distanceLag3 lagdistance
	rename distanceLead1 leadgroup
	rename distanceLead2 leadperiod
	rename distanceLead3 leaddistance

	drop if laggroup == 0
	* Don't save, immediately use 

	********************************************************************************
	*** Find best fit per group

	preserve
	keep laggroup lagperiod lagdistance
	drop if mi(laggroup)
	bysort laggroup: egen minlagdistance = min(lagdistance)
	keep if lagdistance == minlagdistance
	drop minlagdistance
	rename laggroup group
	save group_lag_mindistance_temp.dta, replace

	restore
	keep leadgroup leadperiod leaddistance
	drop if mi(leadgroup)
	bysort leadgroup: egen minleaddistance = min(leaddistance)

	keep if leaddistance == minleaddistance
	drop minleaddistance
	rename leadgroup group
	save group_lead_mindistance_temp.dta, replace

	clear
	use group_lag_mindistance_temp.dta
	merge 1:1 group using group_lead_mindistance_temp.dta

	drop _merge
	bysort group: gen mindistance = min(leaddistance, lagdistance)

	replace lagdistance = . if lagdistance!= mindistance
	replace leaddistance = . if leaddistance!= mindistance

	replace lagperiod = . if lagdistance == .
	replace leadperiod = . if leaddistance == .

	replace lagperiod = -1*lagperiod
	gen period = min(lagperiod, leadperiod)
		* Min because missing is infty 

	keep group period mindistance
	drop if mi(group)
	save groupleadlags_temp.dta, replace

	************************************************************************
	*** Run Synthetic controls on fit lead/lags

	cd /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/botox/synth/100pct/placebos
	use growthgroupspanel_temp.dta, clear 

	*** Merge leads and lags
	merge m:1 group using groupleadlags_temp.dta
	drop _merge mindistance


	** Store group levels in macro
	levelsof group, local(levels)


	*** Replace generate with fitted lead/lags
	tsset group dt
	sort group dt
	gen out = .


	foreach g of local levels{
		disp `g'
		* Treated group gets no time shift
		if `g' == 999{
			replace out = pmt_amt if group == `g'
			continue
			}

		** Grab the stored shifter 
		qui summarize period if group == `g'
		local shift = r(min)

		** If positive, create lead 
		if `shift' > 0{
			replace out = F`shift'.pmt_amt if group == `g'
		}
		** If negative, create lag 
		if `shift' < 0{
			local shiftvalue = abs(`shift')
			replace out = L`shiftvalue'.pmt_amt if group == `g'
		}

		if `shift' == 0{
			replace out = pmt_amt if group == `g'
		}
	}


	** Cut missing data from front of time series
	** Which is limited by the highest lag 
	** Data start in 02
	summarize period
	drop if dt < mofd(mdy(1,1,2002)) + abs(min(r(min), 0))

	*** Cut missing data from end of time series
	*** Which is limited by furthest lead
	*** Data end 9/1/2015 because of second DRG change
	summarize period
	drop if dt > mofd(mdy(9, 1, 2015)) - abs(max(0,r(max)))


	*** Filing date
	disp mofd(mdy(6,5,2007))

	**** Run Synthetic Controls
	synth out out, trunit(999) trperiod(569) figure

	*** Store fitted values
	matrix fit = e(Y_synthetic)
	keep if group == 999
	svmat fit 
	save synth_fit_temp.dta, replace

	********************************************************************************
	*** Compute deterrence
	cd /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/botox/synth/100pct/placebos

	use growthgroupspanel_temp.dta, clear
	keep if group == 999 

	merge 1:1 dt using synth_fit_temp.dta

	local treatment_period = mofd(mdy(6,5,2007))
	disp `treatment_period'

	** Produce difference in series
	gen deterrence = fit1 - out if dt >= `treatment_period'


	** Time discount the difference at 10% 
	local r = 1.1^(1/12)

	gen deterrence_disc = deterrence*1/(`r'^(dt-`treatment_period'))
	replace deterrence_disc = . if dt-`treatment_period' >= 60
	count if !mi(deterrence_disc)
	assert(r(N)==60)

	*** Compute total
	collapse (sum) deterrence_disc
	qui summarize deterrence_disc
	local deterrence = r(min)

	matrix deterrence = (deterrence\ `p', `deterrence')

	rm synth_fit_temp.dta 
	rm groupleadlags_temp.dta 
	rm group_lead_mindistance_temp.dta 
	rm group_lag_mindistance_temp.dta 
	rm growthgroupspanel_temp.dta 

}

clear 
svmat deterrence
rename deterrence1 group 
rename deterrence2 deterrence
sort deterrence
save deterrence_dist.dta, replace


*** Empirical CDF 
cd /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/botox/synth/100pct/placebos
use deterrence_dist.dta 
count 
	* 93
count if deterrence < -41668716
	* 3






